Raport powstał na bazie zbioru ponad 50 tysięcy rekordów z kilkudziesięciu ostatnich lat, dotyczących rozmiarów łowionych w tym okresie śledzi. Dane zostały uzupełnione o brakujące wartości, poddane analizie pod kątem korelacji poszczególnych zmiennych i wykorzystane do zbudowania przez uczenie maszynowe metodą Random Forest regresora pozwalającego przewidywać rozmiar złowionej ryby na podstawie parametrów środowiskowych.
Przeprowadzone badanie wykazało, że istnieje w szczególności silne powiązanie pomiędzy stopniowym wzrostem mierzonych temperatur przy powierzchni wody, a trendem spadkowym w długości łowionych śledzi.
Wczytanie bibliotek dplyr (przetwarzanie danych), ggplot2 i plotly (prezentacja graficzna), easyGgplot2 (rozszerzenia do ggplot2) oraz caret. Dodatkowo wykorzystanie doMC do przetwarzania równoległego.
library(dplyr)
library(plotly) # devtools::install_github(“ropensci/plotly”)
library(ggplot2)
library(easyGgplot2) # devtools::install_github("kassambara/easyGgplot2")
library(reshape2)
library(caret)
library(doMC)
registerDoMC(cores = 4)
raw.data <- read.csv('sledzie.csv', na.strings = '?')
Dane zawierają sporadycznie występujące puste wartości, które zostaną wypełnione medianami występującymi dla odpowiednich wartości kolumn xmonth i recr, które nie zawierają wartości pustych. Uznano to rozwiązanie za akceptowalne dla celów tej pracy, ponieważ liczności zbiorów wartości poszczególnych kolumn są niewielkie (w granicach 50 unikalnych wartości) w porównaniu z liczbą wierszów całej tabeli (ponad 50 tysięcy). Brakujące wartości są na tyle nieliczne, że wpływ ewentualnych przekłamań na wynik pracy można traktować jako pomijalny.
processed.data <-
raw.data %>%
rowwise() %>%
mutate_all(
funs(
replace(., is.na(.), (
function(rownum, v, ref_recr, ref_xmonth) {
if (!is.na(v)) {
return(v);
}
else {
filtered = raw.data %>% filter(recr == ref_recr, xmonth == ref_xmonth)
med = (filtered %>% summarize(median(., na.rm = TRUE)))[[1]]
return(med)
}
}
) (X, ., recr, xmonth)
)
)
)
Zbiór zawiera 52582 rekordów. Poniższa tabelka zawiera podstawowe statystyki dotyczące zbioru danych uzupełnionego o brakujące wartości komórek.
knitr::kable(summary(processed.data[, 1:8])); knitr::kable(summary(processed.data[, 9:16]))
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | |
|---|---|---|---|---|---|---|---|---|
| Min. : 0 | Min. :19.0 | Min. : 0.0000 | Min. : 0.0000 | Min. : 0.000 | Min. : 5.238 | Min. : 0.3074 | Min. : 7.849 | |
| 1st Qu.:13145 | 1st Qu.:24.0 | 1st Qu.: 0.0000 | 1st Qu.: 0.2778 | 1st Qu.: 2.469 | 1st Qu.:13.427 | 1st Qu.: 2.5479 | 1st Qu.:17.808 | |
| Median :26290 | Median :25.5 | Median : 0.1111 | Median : 0.7012 | Median : 5.750 | Median :21.435 | Median : 7.0000 | Median :24.859 | |
| Mean :26290 | Mean :25.3 | Mean : 0.4440 | Mean : 2.0257 | Mean : 9.994 | Mean :21.219 | Mean : 12.7967 | Mean :28.422 | |
| 3rd Qu.:39436 | 3rd Qu.:26.5 | 3rd Qu.: 0.3333 | 3rd Qu.: 1.7936 | 3rd Qu.:11.500 | 3rd Qu.:27.193 | 3rd Qu.: 21.2315 | 3rd Qu.:37.232 | |
| Max. :52581 | Max. :32.5 | Max. :37.6667 | Max. :19.3958 | Max. :75.000 | Max. :57.706 | Max. :115.5833 | Max. :68.736 |
| fbar | recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|---|
| Min. :0.0680 | Min. : 140515 | Min. :0.06833 | Min. : 144137 | Min. :12.77 | Min. :35.40 | Min. : 1.000 | Min. :-4.89000 | |
| 1st Qu.:0.2270 | 1st Qu.: 360061 | 1st Qu.:0.14809 | 1st Qu.: 306068 | 1st Qu.:13.60 | 1st Qu.:35.51 | 1st Qu.: 5.000 | 1st Qu.:-1.89000 | |
| Median :0.3320 | Median : 421391 | Median :0.23191 | Median : 539558 | Median :13.86 | Median :35.51 | Median : 8.000 | Median : 0.20000 | |
| Mean :0.3304 | Mean : 520366 | Mean :0.22981 | Mean : 514973 | Mean :13.87 | Mean :35.51 | Mean : 7.258 | Mean :-0.09236 | |
| 3rd Qu.:0.4560 | 3rd Qu.: 724151 | 3rd Qu.:0.29803 | 3rd Qu.: 730351 | 3rd Qu.:14.16 | 3rd Qu.:35.52 | 3rd Qu.: 9.000 | 3rd Qu.: 1.63000 | |
| Max. :0.8490 | Max. :1565890 | Max. :0.39801 | Max. :1015595 | Max. :14.73 | Max. :35.61 | Max. :12.000 | Max. : 5.08000 |
Wartości średnie oznaczone liniami przerywanymi w kolorze czerwonym.
plot.var <- function(varname, ...) {
ggplot(processed.data, aes(x=get(varname))) +
geom_histogram(color="black", fill="white", ...) +
geom_vline(aes(xintercept=mean(get(varname))), color="red", linetype="dashed", size=1) +
xlab(varname)
}
ggplot2.multiplot(
plot.var("length", binwidth = 0.5), plot.var("cfin1", bins = 30), plot.var("cfin2", bins = 15),
plot.var("chel1", bins = 30), plot.var("chel2", bins = 30), plot.var("lcop1", bins = 30),
plot.var("lcop2", bins = 30), plot.var("fbar", binwidth = 0.08), plot.var("recr", bins = 15),
plot.var("cumf", bins = 9), plot.var("totaln", bins = 12), plot.var("sst", binwidth = 0.15),
plot.var("sal", binwidth = 0.012), plot.var("nao", binwidth = 1),# plot.var("xmonth", binwidth = 1),
cols = 3
)
cor_matrix = round(cor(processed.data), 2)
#cor_matrix[lower.tri(cor_matrix)] <- NA
melted_cor_matrix <- melt(cor_matrix, na.rm = TRUE)
ggplot(data = melted_cor_matrix, aes(Var2, Var1, fill = value)) +
geom_tile(color = "white") +
scale_fill_gradient2(low = "#91bfdb", high = "#fc8d59", mid = "#ffffbf",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Korelacja\nPearsona") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1)) +
coord_fixed()
Mapa korelacji pozwala na wyciągnięcie pewnych wstępnych wniosków dotyczących przydatności poszczególnych zmiennych do budowy regresora przewidującego długość łowionych ryb (atrybut
length).
xmonth możemy w rzeczywistości uznać za atrybut nominalny i sama w sobie nie niesie ona żadnych informacji mogących pomóc w predykcji - jest to atrybut wtórny, będący jedynie niepełną częścią informacji o chronologii.lcop1 i chel1, lcop2 i chel2 oraz cumf i fbar wykazują silną współliniowość. Rozsądne jest pozostawienie tylko po jednym atrybucie z każdej z tych par, ponieważ w przeciwnym razie zaburzone mogłyby być dane o istotności poszczególnych cech.length oznaczającym długość złowionej ryby są atrybuty:
sst - temperatura przy powierzchni wody (korelacja -0.45)nao - oscylacja północnoatlantycka (korelacja -0.26)sst i nao posiadają przy tym wyraźną dodatnią korelację z atrybutem X, co zgodnie z założeniem o zachowaniu przez dane wejściowe chronologii mogłoby sugerować, że potencjalnie znaleziony regresor może uznać te zmienne za istotne czynniki predykcji.processed.data <- processed.data %>% select(-lcop1, -lcop2, -cumf, -xmonth)
Poniższy wykres przedstawia zależność długości ryby (mierzonej z dokładnością ) od czasu (numeru sekwencyjnego próbki w zbiorze). Dodatkowo, próbki zostały pokolorowane w zależności od sst (temperatury przy powierzchni wody) - jaśniejszy odcień koloru niebieskiego oznacza wyższą temperaturę, a u dołu wykresu znajduje sie pasek pozwalający na przełączanie pomiędzy próbkami w momentach, gdy oscylacja północnoatlantycka (atrybut nao) osiągał wartości dodatnie, a tymi wykonanymi przy ujemnej wartości tej zmiennej (nao_range odpowiednio 1 i -1).
Na wykresie ze względu na wydajność nakreslono co 35 rekord ze zbioru.
ggplotly(
ggplot(processed.data[seq(1, nrow(processed.data), freq), ] %>% mutate(nao_range = sign(nao)), aes(X, length)) +
geom_point(aes(color=sst, frame=nao_range)) +
labs(x = "numer sekwencyjny próbki",
y = "długość ryby (cm)",
color = "Temp. pow.\nwody (st. C)",
title = "Długość ryby, czas i temperatura")
)
Wizualizacja zdaje się być kolejnym czynnikiem mogącym potwierdzać przypuszczenie o związku temperatury przy powierzchni wody z długością łowionych ryb. W początkowym okresie badania (około 17000 próbek) wartość temperatury wykazywała tendencję spadkową (zilustrowane ciemniejącym odcieniem błękitu punktów), rosły zaś wskazania pomiaru długości ryb. W dalszej części wykresu temperatura przy powierzchni wody rośnie, spadają zaś pomiary długości śledzi.
Trudniejszy do oszacowania wydaje się być wpływ parametru oscylacji północnoatlantyckiej. Jest to zjawisko cykliczne, która to cykliczność polega na okresowej dominacji odczytów dodatnich lub ujemnych - co zresztą znajduje potwierdzenie w wykresie, jeśli spojrzeć na okresowe zagęszczenie próbek dla ujemnego i dodatniego nao. Atrybut ten jest też w dość dużym stopniu skorelowany z sst - w okresie dominacji ujemnych odczytów nao dominują też niższe wartości sst i większe długości łowionych śledzi.
lengthPodział zbioru w proporcjach: 60% zbiór uczący, 20% zbiór testowy, 20% zbiór walidacyjny.
set.seed(23)
inTraining <- createDataPartition(y = processed.data$length, p = 0.6, list = FALSE)
training <- processed.data[inTraining, ]
testing.and.validation <- processed.data[-inTraining, ]
inTesting <- createDataPartition(y = testing.and.validation$length, p = 0.5, list = FALSE)
testing = testing.and.validation[inTesting, ]
validation = testing.and.validation[-inTesting, ]
Uczenie z wykorzystaniem wielokrotnej oceny krzyżowej.
ctrl <- trainControl(
method = "repeatedcv",
number = 2,
repeats = 5
)
set.seed(23)
fit <- train(length ~ ., data = training, method = "rf", trControl = ctrl, verbose = TRUE, importance = TRUE)
## Loading required package: randomForest
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
##
## margin
## The following object is masked from 'package:dplyr':
##
## combine
pred <- predict(fit, validation)
measures = postResample(pred = pred, obs = validation$length)
rmse <- measures[[1]]
rsquared <- measures[[2]]
Wygenerowany regresor dla zbioru walidującego osiągnął RMSE=1.1204929 i R2=0.5384306. Biorąc pod uwagę zakres wartości przewidywanej zmiennej length, taką wartość pierwiastka błędu średniokwadratowego można uznać za zadowalającą przynajmniej jako początkowy punkt odniesienia do analizy danych. W przypadku osiągniętej wartości R2 wyliczony model wyjaśnia około 50% wariancji zmiennej length i około 30% jej odchylenia standardowego. Nie jest to zatem predyktor zapewniający bardzo wysoką pewność otrzymania prawidłowego rezultatu, niemniej jednak warto wyciągnąć także wniosek z poniższego wykresu.
ggplot(validation %>% mutate(predicted_length = pred), aes(X, length)) +
geom_point(aes(color = predicted_length - length)) +
scale_colour_gradientn(colours = c("red","green", "blue")) +
labs(x = "numer sekwencyjny próbki",
y = "długość ryby (cm)",
color = "Błąd predykcji",
title = "Relacja długości ryby do wartości przewidzianej przez model")
Wykres ten ponownie przedstawia zmienne w czasie zmierzone długości ryb, jednak tym razem punkty zostały pokolorowane stosownie do błędu predykcji dla danej próbki. Uzyskany model jest dosyć konserwatywny i niechętnie wychodzi poza pewien zakres wartości (na wykresie kolor zielony), notując w nim wysoką dokładność, zaś myląc się głównie w przypadkach skrajnych. Można przypuszczać, że liczność tych ostatnich jest główną przyczyną dosyć przeciętnego rezultatu współczynnika determinacji.
importance <- varImp(fit)
## Loading required package: randomForest
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
##
## margin
## The following object is masked from 'package:dplyr':
##
## combine
plot(importance)
Analiza ważności poszczególnych zmiennych potwierdza wykorzystanie przez predyktor w największym stopniu zmiennej
sst; druga co do bezwzględnej wartości korelacji zmienna nao nie została zaś w ogóle wykorzystana. Pozwala to na potwierdzenie przyjmowanego dotychczas przypuszczenia, że wzrost temperatury przy powierzchni wody ma kluczowe znaczenie dla obserwowanego trendu malejących wymiarów łowionych śledzi.